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ABSTRACT 

. We propose a method for determining the two-dimensional orbits of interstellar matter (ISM) 

I in galactic disks using observed line-of-sight velocity fields. Although four components, consisting 

of two of position and two of velocity, are necessary to determine two-dimensional orbits, the 
^ , velocity perpendicular to the line of sight can hardly be determined from observations. We 

I show that the ISM orbits can be approximately determined by assuming that (1) the angular 

momentum is conserved along each orbit, (2) all orbits are closed in a single rotating frame, and 
(3) the system is symmetric through the origin. Since the angular momentum can be estimated 
on the major axis of a disk, a velocity perpendicular to the line of sight can be determined 
^ ■ according to the conservation of angular momentum, and thus two components of velocity are 

'sj" \ obtained. The ISM orbits can be calculated by solving differential equations given by an obtained 

' two-dimensional velocity field when a certain pattern speed is given. The typical deviation from 

. the correct orbit is less than the observational resolution according to our test calculations for 

\ modeled velocity fields. We applied this method to the observed data of NGC 4569 and obtained 

^'j^ ' possible orbits for this galaxy. Furthermore, we calculated the mass distribution using our result. 
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^ ' 1. Introduction 



Although orbits of stars or interstellar matter (ISM) on galactic planes are fundamental in studying 
galactic dynamics, accurate shapes of these in real galaxies are difRcult to discern because two-dimensional 
\ motion can hardly be observed. Observations of line emission give only the line-of-sight velocities of the 

ISM in spiral galaxies. Hence, we often discuss galactic dynamics assuming that the ISM motion is circular, 
since the velocity fields of galaxies often show so-called spider diagrams, which can be reproduced by circular 
rotation around their galactic centers. Obtained rotation curves provide mass distributions in spiral galaxies 
(Sofue et al. 1999, 2003a). Barred galaxies, however, show strong noncircular motions, and approximations 
of circular orbits fail. For these galaxies, the computer simulation of ISM dynamics has played an important 
role in studying galactic dynamics, and many models have reproduced observations fairly well (Huntley 1978; 
Wada et al. 1998). However, the link between simulations and observations is generally difRcult, since there 
are many unidentified parameters, such as rotation curve, bar strength, pattern speed, viscosity, etc. And 
although the shape of a galactic bar is often expressed in the form $h(i?) cos (mcf)), where R, to, and <j> denote 
galactic radius, integer, and the angle around the galactic center, respectively, it has not been determined 
by observations. 

Therefore, it is desirable to determine two-dimensional ISM motions on the galactic plane directly from 
observations in order to model real galaxies. While Kuno & Nakai (1997) inferred an ISM orbit by tracing 
the velocity vectors at the spiral phase using the CO data of M51 and found that the orbit deviates greatly 
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from a purely circular motion, there are few galaxies for which two-dimensional orbits have been derived 
from observational data. 

In order to trace two-dimensional motion in a galactic disk wc need four values, two of which arc of 
position and two of which are of velocity. Both position components are directly given by observation. 
A velocity component perpendicular to the line of sight cannot be determined, although the line-of-sight 
component is given by observations. 

In this paper, we propose a method for determining the two-dimensional orbits of ISM on galactic planes 
directly from observations on the basis of simple assumptions. Wc applied this method to model velocity 
fields to confirm that it works well. We then applied this method to a real galaxy, NGC 4569, for which 
^^CO(J =1 — 0) data has been obtained with Nobeyama Millimeter Array (NMA) observations. 



When the potential of the galactic disk is considered to be dominated by an axisymmctric component 
and bar or spiral components are small enough, the force acting on the ISM is directed approximately 
centripetally and the angular momentum is nearly constant. In this section, we describe a method for 
determining an ISM orbit on a galactic plane based on the assumptions that (1) the angular momentum is 
conserved along an orbit, (2) all orbits are closed in a single rotating frame, and (3) the system is symmetric 
through the origin. 

We choose a Cartesian coordinate whose {x, i/)-plane coincides with a galactic plane. In viewing the 
galaxy face-on, we let the x- and y-axes coincide with the major and minor axes of the disk, respectively. 
The galactic center is located at {x,y) — (0,0). The velocity components parallel to the x- and y-axes are 
denoted by u and v, respectively. When the galactic disk is inclined by i, the projected coordinate is defined 
by Xproj, yproj, Mproj, f pro j , whcrc Xproj = X, i/proi = V^osi, Uproj = u, Vproj = vsini. The obscrvablc hne-of- 
sight velocity corresponds to Vproj- While three components {x,y,v) are known from observation with the 
inclination i, one component, u, is unknown. 

Let a galaxy rotate counterclockwise and the lower side of the disk {y < 0) be near us. First, we choose 

two points on the major axis so that they are symmetric with respect to the galactic center; these points are 
defined as (a;J,2/o) = (<i)0) and (x^jy^) = (— a,0). The superscript denotes a quadrant of position; 1, 2, 3, 
and 4 mean the first, second, third, and fourth quadrants, respectively. The subscript represents position; 
denotes an initial point and 1 denotes a neighboring point. Here we mainly think of two trajectories in the 
first and the second quadrants, that is, y > 0. We omit a superscript when we do not specify the quadrant. 

The velocities Vq and Vq at these points are determined from observation. We set Vq = —Vq = vq by 
averaging the observational data at two points in order to reduce observational errors. We suppose that one 
orbit crosses these two points and that the angular momentum, L, is constant along this orbit. Since the 
angular momentum in the Cartesian coordinate is 



and y = along the major axis, the angular momentum along this orbit is determined uniquely as i = avQ. 
Since v at each point (a;, y) is given by observation, u can be estimated to be 



2. 



Method 



L = XV — yu 



(1) 



u = {xv — L)/y 



(2) 



at each point. 
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Second, wc scit the; velocity uq to be wo/tan^ at {xq, yo) using a free parameter, 6. Then a neighboring 
point {xi,yi) is determined in the vector form approximately by 

Xi = xo + (vo - Vpat)di, (3) 

where X and v denote position and velocity vectors, respectively The vector Vpat is the pattern speed, whose 
absolute value is proportional to the radius and whose direction is always tangential to the position vector. 
The constant dt is the time interval. This equation is a formula of the Euler method for solving a differential 
equation numerically. The time interval dt is chosen so that each subsequent point is separated from the 
initial point by half the spatial resolution of the observation, based on the Nyquist sampling. We allow dt 
to take a negative value; positive and negative dt correspond to trajectories traveling counterclockwise and 
clockwise from the initial points. Two trajectories, therefore, start from one initial point. 

The vector component Vi at the next position (xi, j/i) can be determined by observation, and the other 
component, ui, may be estimated by the conservation of angular momentum. We calculate an orbit in 
the first quadrant simultaneously with the third quadrant, which is the same as in the second and fourth 
quadrants. Then, to reduce the observational noises, we average velocity components v at symmetric points 
relative to (0,0) every calculation step; for example, v} = (v} — vf)/2. If data is lacking at the next point, 
wc dcifine the value at this point to be the same as that of the previous; Vfe = Vfc_i, where k is any integer. 
We repeat this process until an orbit reaches x = 0. For the next stage, we vary a parameter 9 by A9, which 
was taken to be 1° in this paper, so that orbits from (±a, 0) meet one another at x = 0. We avoid a solution 
whereby a radius of the orbit becomes extremely large, |xfe| > Tupper- In this paper we adopted rupper = Sxq. 

Figure 1 presents a flow chart and a scheme for this method. Although we explain our method by using 
the Euler method for simplicity, we actually use a fourth-order Runge-Kutta method to calculate orbits more 
accurately. 

— Figure 1 — 

3. Application To Observed Data 
3.1. Data 

In this paper we applied this method to NGC 4569 ^^C0( J =1 — 0) data, which was taken from the 
high-resolution Virgo CO survey (Sofue et al. 2003b). The adopted parameters and observational parameters 
are listed in Table 1. The inclination and position angles were determined by applying an elliptical fitting 
to a K-band image (Jarrett et al. 2003) using the task ellipse in the IRAF STSDAS software. The position 
angle steeply decreases toward the center within R < 3'.'5, which reminds us of a central bar. The position 
angle largely increases outward beyond R > 131", which reminds us of a warping of the stellar disk. In order 
to avoid these effects, we estimated inclination and position angles by averaging them between R = 3'.'5 
and 131". Consequently, we obtained the position angle 18° and the axial ratio q = 0.52, which yields an 
inclination of i = 63°;, according to the equation of TuUy (1988), 

i = 3°+ cos" V(3^ - 0.22)/(l - 0.22). (4) 

NGC 4569 shows an S-shaped velocity field, which indicates the clear existence of noncircular motion because 
the isovelocity contour of the systemic velocity is not parallel to the minor axis at all. The Digitized Sky 
Survey (DSS) B-band and CO images of this galaxy are shown in Figures 2a and 2b. The contours in Figure 
2b represent the velocity field. 
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According to Tremaine & Weinberg (1984), the pattern speed can be estimated by using the intensity- 
weighted line-of-sight velocity and the position of the major axis based on the equation of continuity. We 
apphed this method to this galaxy using the /T-band image and the CO velocity field. The intensity at 
the center is quite large, which reminds us of nuclear activity and does not seem to satisfy the equation of 
continuity. Therefore, we estimated the pattern speed in the range 2" < R < 9" and obtained a pattern 
speed of 74 km s~^ kpc~^. Note that this value is tentative, since gaseous velocity does not always coincide 
with the stellar velocity. 

— Table 1 — 

— Figure 2 — 

3.2. Results 

We applied our method to this galaxy to determine orbits on the galactic planes. For NGC 4569, we 
fixed the step of calculation at 1'.'55, which is half of the beam size. We calculated less than 12 orbits that 
cross the points {x,y) = (±0'.'775n, 0), where n is an integer, on the major axis. We show the resultant 
orbits viewed face-on in Figure 2c. We imposed these orbits on the CO imago, which is shown in Figure 
2d. In NGC 4569, orbits tends to be strongly elongated in the direction of the major axis. The axes of the 
orbits are slightly inclined in the counterclockwise direction from the major axis. This large deviation from 
circular orbits causes the strong disturbance in the velocity field. 

Since the estimated pattern speed is tentative, we calculated orbits taking other pattern speeds (34, 54, 
and 94 km s~^ kpc~^) to find how the resultant orbits depend on the pattern speed. Results are shown in 
Figure 3. Fewer orbits can be obtained when a larger pattern speed is adopted. In any case, resultant orbits 
are strongly elongated in the direction of the major axis and are inclined in the clockwise direction from the 
major axis. 

— Figure 3 — 

Since characterizing the presence of noncircular motions is sensitive to the values of inclination and 
position angle adopted (Schoenmakers et al. 1997), we calculated orbits taking other inclinations and 
position angles that differ by 3° from i = 63° and P. A. = 18°. Results are shown in Figure 4. Only a few 
orbits can be obtained for some adopted inclinations and position angles. In these cases also, resultant orbits 
are elongated and are inclined in the counterclockwise direction. 

— Figure 4 — 

While some orbits seem to cross others slightly, the gas orbits should not cross each other, otherwise the 
gasses would collide, shock, and lose angular momentum. Since the amplitude of the orbits' intersections is 
less than beam size, these intersections are thought to be due to a calculation error. 
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4. Error Estimation 

Below we have applied our method to model velocity fields in order to confirm that it works properly 
and to estimate errors. Consequently, we conclude that our calculation works well within a beam size and 
that the calculation gives the best solution when we adopt a half beam size or a beam size as a step size. 

4.1. Model 1: Angular Momentum Is Conservative And Orbits Are Aligned 

First, we made a simple model in which each orbit is an ellipse whose axial ratio is 2 and the position 
angle of each orbit is the same as the others. The velocity at the major axis is fixed at 200 km s~^, and the 
angular momentum is thoTight to be constant along each orbit. The inclination of the disk is fixed at 60°. 
The velocity field is given in an infinite resolution. When step sizes of 2". 88 and 0".288 are adopted, typical 
deviations from the correct orbits are less than 1 /20 of a step size. We conclude that the calculation works 
well for this ideal case. 



4.2. Model 2: Angular Momentum Is Conservative And Orbits Are Misaligned 

Next, for a more complicated case, we made a model in which each orbit is an ellipse whose axial ratio 
is 2 and the position angle of each orbit increases with radius; the semimajor radius is 0."288n (n being an 
integer), and the position angle increases by 0.°5 when the semimajor radius increases by 0."288. Angular 
momentum is constant along each orbit. The inclination of the disk is fixed at 60°. Velocity fields are given 
by smoothing with a Gaussian whose full width at half-maximum (FWHM) is 2. "88 on the basis of a typical 
beam size of the NMA Virgo CO survey. We calculated orbits with various step sizes. An averaged deviation 
between 2. "88 and 10" gets the smallest value, 0."9, when we adopt a beam size as a step size. 

4.3. Model 3: Angulcir Momentum Is Not Conservative And Orbits Are Misaligned 

It is not precise to assume that angular momentum is constant along an orbit in a field in which the 
gravitational force is not always centripetal. To see the effect of the variation of angular momentum along 
an orbit, we applied our method to a model of orbits in a weak bar potential. We took a model from Wada 
(1994). We chose a bar strength e = 0.05 and a damping rate A = A/kq = 0.1. The angular momentum 
varies by ~ 10%. We assume that the rotational velocity reaches a maximum of 203 km s~^ at a radius of 
1.6 kpc. The pattern speed, Opat, is fixed at 70 km s~^ kpc~^. The distance is taken to be 16.1 Mpc, which 
is the distance to the Virgo Cluster. The averaged deviation is minimized at 0."5 when we adopt a step size 
equal to a half of the beam size. 

If the angular momentum varies by a%, a velocity vector has an error of a% at maximum. When we 
fixed the step size at half a beam size, the position error of the orbit is a% of half a beam size at maximum 
for each step. An error due to the assumption of a constant angular momentum would not be negligible 
when the step number is large. 
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4.4. Step Size And Error 

We found that deviatfon is the smallest when the step size is fixed at a beam size or a half beam size. 
This is interpreted as follows: The deviation grows as the step number increases. On the other hand, the 
smaller the step size is, the more accurately the Rungc-Kutta integrator works. Since the observed data are 
smeared with the beam size, we cannot get information for a scale smaller than the beam size. When the 
step size is smaller than the beam size, the calculation error only increases. Therefore, the step number and 
step size should be as small as possible and the step size should not be much smaller than the beam size in 
order to reduce the error. Consequently, the deviation is the smallest when we adopt a beam size or half a 
beam size as the step size. 

4.5. Dependency On The Pattern Speed 

In order to investigate the dependency on the adopted pattern speed, we applied our method to model 
3, assuming incorrect pattern speeds of 30, 50, and 90 km s^^ kpc^^. The averaged deviations are l."5, 
0."7, and 0."9, respectively. The deviation increases for an incorrect pattern speed. The averaged deviation 
is, however, smaller than the beam size for this model. 

5. Mass Distribution 

One of the most useful applications of our method is to determine the mass distribution in a spiral galaxy 
even if the rotation curve cannot be obtained because of strong noncircular motion. Following Binney & 
Tremaine (1987), we introduce a frame rotating at ilp^t = (0, 0,npat) and effective potential 

$eff = $-^^^^at^^ (5) 

where $ is a potential in a rest frame and i? is a radius. In this frame, the equation of motion becomes 

i^ = V$cfr-2(f2pat xr). (6) 

Here the acceleration r is known from the resultant orbits and two-dimensional velocity vectors. A mass 
density can be calculated according to the Poisson equation 

V^$(r) = 47rGp(r). (7) 

We found that the gravitational vector is oriented to the galactic center within 10° by calculating it 
from the obtained orbits. Therefore, the potential is nearly spherical and is expressed as 

a2$(r) /dx'^ ~ 52$(r)/a2/2 _ d^^{r)/dz'^ (8) 

with a function of only radius r. Therefore, we assumed 

d'^^r)/dz'^ ~ {d'^^r)/dx'^ + a2$(r)/V)/2 (9) 

in order to calculate the mass distribution in NGC 4569 using our result. We calculated the mass distribution 
with pattern speeds of 34, 54, 74, and 94 km s~^ kpc~^. 
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FigTirc 5 shows the mass distribution on the galactic plane of NGC 4569. In this paper, wc calculated 
mass densities inside the largest obtained orbits. NGC 4569 shows an elongated distribution which is inclined 
in the counterclockwise direction from the major axis. The mass distribution near the borderlines is disturbed 
because there are not many data points to calculate. 

A typical error arising from a 10.4 km s~^ error in the line-of-sight velocity is 5.7 Mq pc~^ when we 
suppose (1) a velocity vector to be parallel with y-axis; (2) an error in a gradient of the potential to be 
(5|V<i>| ~ 2(5w/(5t + 2ripat<5w ^ 2v5v/5y + 2Vlp^^5v, and (3) an error in mass density to be AirGSp ^ 25\V^\/ 5y. 
Here, we adopt 6v = 10.4 km s^^, fipat = 74 km s^^ kpc^^, v = 100 km s^^, and 6y = 1.55 x 78 pc. 

— Figure 5 — 

6. Remarks In Applying Our Method 

We presented the ISM orbits as being approximately determined in assuming that (1) the angular 
momentimi is conserved along each orbit, (2) all orbits are closed in a single rotating frame, and (3) the 
system is symmetric through the origin. Our method is useful for modeling a real galaxy because it is 
calculated directly from observational data. However, we must note the following points: First, our method 
is plausible for a potential that is dominated by an axisymmetric component and in which the loss of angular 
momentum due to the shock of the gas is negligible. Otherwise, the angular momentum is not conserved 
and our method fails. Second, this method cannot be applied to a patchy velocity field. Third, this method 
cannot work at the corotation radius where the rotational velocity is equal to the pattern speed. Fourth, 
resultant orbits can deviate by observational resolution. 
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Table 1: Adopted Parameters and Observational Parameter 



Name 
R.A.(J2000) 
Dec.(J2000) 
Morphology type 
inclination 
position angle 
systemic velocity [km s^-'^] 
beam size 
r.m.s. noise [mjy] 
velocity resolution [km s^^] 



NGC 4569 
12^36'"49.80''*i 
+13°09'46.3"*i 
SAB(rs)ab*2 
63°*3 

-235 ± 4*2 
4."5 X 3."r4 
25*4 
10.4*4 



*i The 2MASS Large Galaxy Atlas (Jarrett et al. 2003) 
*2 RC3 (de Vaucouleurs et al. 1992) 

*^this work (elliptical fitting using the task ellipse in the IRAF STSDAS software for 2MASS K-band image) 
*4Sofue et al. (2003a) 
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C Xq. Vo, Vq ^ 



Xq. Vo . Uo, Vo 



e = e+A e 




repeat (*) 



c 




(b) 



A y-axis 
minor axis 




(xio,yio)=(a, 0) 



average velocities 



Fig. 1. — A flow chart (left) and a scheme (right) for our method. 
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Fig. 2.— (a) DSS B-band image of NGC 4569. A square indicates the area of (b). (b) A ^^CO (J=l-0) map 
with isovelocity contours; -335, -285, -235, -185, and -135km s^^. (c) Resultant orbits and two-dimensional 
velocity vectors viewed face-on. (d) Resultant orbits superposed onto the CO image of (b). 
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Fig. 3. — Resultant orbits and two-dimensional velocity vectors assuming that the pattern speed is 34, 54, 
74, and 94 km s~^ kpc~^ 
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Fig. 4. — Resultant orbits and two-dimensional velocity vectors assuming that the inclination is 60°, 63°, 
and 66°, and that the position angle is 15°, 18°, and 21°. 
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Fig. 5.— Mass distribution for NGC 4569. The contour levels are 26, 18, 13, 9.1, 6.4, 4.5, 3.2, 2.3, and, 1.6 
M0 pc-3. 



